Analytical investigation of modulated spin torque oscillators in the framework of 
coupled differential equations with variable coefficients 
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Modulation of Spin Torque Oscillators (STOs) is investigated by analytically solving the time- 
dependent coupled equations of an auto-oscillator. A Fourier series solution is proposed, leading to 
the coefficients being determined with a linear set of equations, from which a Nonlinear Amplitude 
and Frequency Modulation (NFAM) scheme is obtained. In this framework, the NFAM features 
are related to the intrinsic STO parameters, revealing a frequency-dependence of the harmonic- 
dependent modulation index that allows a modulation bandwidth to be defined for these devices. 
The presented results expose a rich parameter space, where the modulation and the STO's operation 
conditions define the observed modulation features. The Fourier-series representation of the time 
signal is suitable for studying periodic perturbations on the auto-oscillator equation. 



I. INTRODUCTION 

Spin torque oscillators^ (STOs) are spin torque^"— 
driven devices whose high-frequency tunability and nano- 
sized dimensions are promising for diverse technological 
applications. The oscillators' functionality can be very 
broad, depending on the specific application and its op- 
eration frequency. Two phenomena are of particular use, 
namely (i) the ability to phase lock or synchronize to 
another periodic signal, and (ii) the possibility of mod- 
ulating a base-band signal.— In the case of STOs, both 
phenomena have been experimentally observed. 

Phase locking of STOs has been extensively studied in 
the past few years. Injection locking of an STO to an 
external source is already well-established in both exper- 
imental^"— and numerical^"— works. Phase locking can 
also be achieved by mutually coupling several STOs. In 
this case, a distinction must be made between nanocon- 
tacti^ and nanopillar— geometries. In nanocontacts, syn- 
chronization mediated by propagating spin waves^^ii^ has 
been satisfactorily achieved both experimentallyi^Ti^ and 
numerically^ for GMR-based STOs. In contrast, the mu- 
tual locking of nanopillars has been mainly successful in 
numerical studies i^"— 

Modulation of STOs has also been shown to be exper- 
imentally straightforward giving rise to several applica- 
tions. STOs nanopillars have been proposed as hard-drive 
read-heads where the stray field from the perpendicular 
media acts a modulating source^^ both in and 
TMR^^ spin valves. From the perspective of communica- 
tion applcations. Frequency Modulation (FM) has been 
observed in GMR nanocontacts^"— and more recently 
in nano-oxide layer (NOL) nanocontact STOs<^ Digital 
communication schemes, in particular Frequency Shift 
Keying (FSK), has been demonstrated in vortex-based 
nanocontact STOs up to the limit of analog FM<^i2i 
Similar modulation effects are most likely to be observed 
in other STO geometries that take advantage of Per- 
pendicular Magnetic Anisotropy (PMA) materials^ and 
nanopillars located in a microstrip resonatori^ 



In the pioneering frequency modulation experiment;^ 
the authors took advantage of the strong phase-power 
coupling of GMR-based STOs to obtain a frequency- 
modulated voltage output from a base-band current tone. 
Although the gross features of frequency modulation were 
present, significant discrepancies from the model were 
observed: (i) an unexpected carrier frequency shift as 
a function of modulation strength, and (ii) asymmet- 
ric sideband amplitudes in contrast to the expected n-th 
order Bessel functions. A mathematical model was pro- 
posed by Consolo et al.— to fit these features. Here, the 
temporal signal s{t) was defined to be both amplitude- 
modulated and frequency-modulated 



s{t) = ^[a;(t)] cos(cl)c + uJi[x{t)])t, 



(1) 



where Uc is the carrier frequency, and the amplitude 
A and instantaneous frequency Ui are polynomial ex- 
pansions of a base-band message x{t). This combined 
action of amplitude and frequency modulation was 
named "Nonlinear Amplitude and Frequency Modula- 
tion" (NFAM). Using this approximation, it was possible 
with very good accuracy to fit the GMR-STO modula- 
tion data with a polynomial expansion up to the third 
order.— However, the information acquired from these 
coefficients does not reflect the intrinsic mechanism of 
NFAM in STOs, nor its consequences for a specific choice 
of experimental conditions and device characteristics. To 
investigate the origin of NFAM in STOs, we here derive 
the modulation spectrum as a function of the intrinsic 
STO parameters. 

The modulation spectrum of an STO is obtained by 
perturbing the auto-oscillator general equatioi*^ with a 
slow time- varying tone. Such perturbation creates cou- 
pled phase and amplitude variations, which lead to 
NFAM. It is shown that a Fourier series gives an ade- 
quate representation of s{t), in which case the problem 
can be reduced to determining the coefficients from a 
linear set of equations. Moreover, the Fourier coefficients 
obtained by this method give quantitative information, 
such as the carrier frequency shift, the modulation index, 
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FIG. 1. Nanocontact pseudo spin valve consisting of a thick 
(fixed) and a thin (free) magnetic layer (gray) decoupled by 
a thin metallic spacer (black). A nanocontact is patterned 
on top of the free layer, allowing the device to be driven 
with a spatially confined current. It is assumed that the auto- 
oscillator Eq. ((2)) is a good approximation for such a structure. 
The free- layer magnetization m is free to precess about the ap- 
plied field Happ- The fixed layer is assumed to be magnetized 
predominantly in plane, and its magnetization Ps subtends 
an angle 70 relative to the free layer magnetization. A steady 
precession is achieved when a spin polarized current exerts 
sufficient spin torque on the free layer to compensate for the 
magnetic damping. 



and the modulation bandwidth of the STO. 

This paper is divided as follows: in Section II, the gen- 
eral formulation of the problem and its solution is given. 
It is shown that a carrier frequency shift appears as a 
consequence of power and phase coupling. Moreover, the 
modulation index is obtained by calculating the power 
spectral density, and shows an unexpected modulation- 
frequency dependence. In Section III, the analytical solu- 
tion is evaluated for different modulation and operation 
conditions, and shows a qualitative agreement with ex- 
perimental observations. Concluding remarks are given 
in Section IV. 



II. MODULATION OF A NONLINEAR 
AUTO-OSCILLATOR 

In this section, the power spectral density (PSD) of a 
frequency-modulated nonlinear auto-oscillator is analyt- 
ically calculated, using the general model proposed by 
Slavin and Tiberkevich;^ 



dc 
di 



+ iuj{p)c + r+(p)c — r_(p)c = 0, 



(2) 



where p is the oscillation's power, c = y/pe^'-'^ is the 
oscillation's amplitude, (f> is its time-dependent phase, 
and aj(p), r+(p), and T^{p) are respectively the power- 
dependent oscillation frequency, damping, and negative 
damping. The power dependencies of these quantities 
are treated in Ref. |4ll as polynomial expansions. We as- 
sume that the model of Eq. ([2]) is a good description of 
the nanocontact geometry schematically shown in Fig. [T] 



Here, a metallic nanocontact is patterned on top of an 
extended GMR spin valve. A strong magnetic field Happ 
is applied to ensure the saturation of the magnetically 
active or "free" layer material. When a bias current is 
driven through the nanocontact, spin torque is exerted 
on the underlying free-layer area, counteracting the mag- 
netic damping action, and leading to a small-angle pre- 
cession about the applied field. Depending on the applied 
field angle, the free layer might exhibit propagating spin 
waves or a so-called localized bullet.—"— In this paper, 
we restrict our attention to a perpendicular applied field 
without loss of generality. 

In order to solve the modulation problem, we recall 
the experiment of Ref. [2^ where an ac current is added 
to the dc bias current. Consequently, we define the slow 
time- varying current as 



Ht) = Idc (1 + /iCOSWm^) , 



(3) 



where /i is the modulation strength, defined as the ratio 
between the ac amplitude and the dc current, and LOm is 
the modulation frequency. We assume that uj,„ ^ '^sto 
(where ujsto is the STO's free-running frequency) in or- 
der to be considered a modulating frequency. On the 
other hand, if tOm ~ ^sto, the STO might instead be- 
come injection-locked. The current is thus included in 
the negative damping parameter, r_(p), which can be 
generally expanded polynomial in power 



r_ (p) = al{t) {qo + qip + q2p' + ...) 



(4) 



where a — eh'y /2eMsV , e the spin-polarization effi- 
ciency, 7 the gyromagnetic ratio, e the electron charge, 
Ms the free layer's saturation magnetization, V the free 
layer's current-carrying volume, and the coefficients qt, 
1 = 0,1, 2, ... are assumed to describe the sample's char- 
acteristics. It is further assumed that lit) creates a power 
perturbation of the form p = po(l + 2(5p), where po is the 
STO's free-running power. This approximation is valid as 
long as the modulation strength is small, which we con- 
sider a common scenario in STO modulation experiments 
(for instance, /.t < 0.1 in both Ref. [2^ and Isol) . 

Separating Eq. ([2]) into real and imaginary parts, and 
expanding in power to first order, we obtain a set of dif- 
ferential equations with variable (time-dependent) coef- 
ficients for the STO's power and phase: 



dSp 
~dt 



^Ci cos U!mt + (^6*2 COS U!mt — Tp) 2Sp, (5a) 



— = UJSTO + 2vTp5p, 
where we define the constants 



Ci = r_(po), 
C2 = ('r_(po) 



dT-{p) 



dp 



Po 



Po 



(5b) 

(6a) 
(6b) 



Here, ujsto — ^o+i^^p, — Jf^o{H— Ms) is the FMR 
frequency, /ig is the vacuum permeability, Tp = aa;o(^— 1) 
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is the restoration rate, a is the Gilbert damping param- 
eter, ^ = Idc/hh is the supercriticality parameter, and 
Ith is the threshold current for spin-torque-driven oscilla- 
tions. The restoration rate is assumed to be constant with 
respect to the free-running oscillations. Qualitatively, it 
gives a measure of how fast the STO reacts to a pertur- 
bation, and so plays a fundamental role determining the 
synchronization speed for these devicesJ^i^i^ 

We propose a Fourier series as a solution of Eq. ([5]) a, 
based on the periodicity of the time-dependent term: 

oo 

6p = Aq + An sin nUJrnt + Bn COS nUJmt. (7) 
n=l 

Introducing this solution into Eq. ([5]) a gives an infinite 
system of equations for the Fourier coefficients. It is pos- 
sible to further reduce the problem to the determination 
of An , leading to a system of linear equations which can 
be solved numerically (see Appendix A). The solution, 
expressed as a single sinusoid, is 

oo 

Sp ^ Aq + ^ y^Bl + Al cos{nuj,nt - V'n), (8a) 

n=l 

(j) = [uJsTO + IvTpAoY 
+ y -^^Bl + Al sin{numt - ^„). (8b) 

, nUJm 
n—1 

The notation of Eq. (jS]) explicitly shows a harmonic- 
dependent phase ipn = arctan(A„/i?„). This phase plays 
a fundamental role in the form of the PSD, as discussed 
below. The first indication of NFAM is visible in the first 
right-hand term of Eq. (j8b[) . where the oscillation fre- 
quency is shifted by 27r/s — 2i'TpAo. From the exact 
solution given in Eq. (jASp . we obtain 



2TTfs = - 



^wC2Bi 



(9) 



Consequently, the shift is directly proportional to the 
modulation strength and to the STO nonlinearities, via 
the constant C2. 

In order to identify the next feature of NFAM — namely 
the asymmetry of the sideband — the PSD must be cal- 
culated. By expressing the sinusoidal functions as expo- 
nentials and performing Taylor expansion (Appendix B), 
the PSD can be expressed as a series of convolutions 

PSD = po5{uJsTO' ) * f(l + ^0)^5(0) 



X X 

+ ^ -^5{nu)m.) + -^5{-nuj^n) 



n=l 

00 



* Y[Mf3n)S{0)+J2 

n=l j=l 

+ {-iyXidi~njLU,n))\ 



\Xn\l 



{Xi5{nju,n) 



(10) 



complex conjugate, and the shifted carrier frequency is 
written as losto' = ^sto + 27r/s. The notation S{xo) = 
6{x — xq) is used, where 6 is the Kronecker delta. The 
product symbol represents in this case a series of convo- 
lutions. The real harmonic-dependent modulation index 
is defined as 



_ 2iyTp\Xn 

Pn — 



where we define the complex variable X„ = i?„ -I- iAn 
(hence V'n = a'rg(X„)), the notation (•) represents the 



(11) 



Comparison with the FM modulation index, h = 
A///m, suggests that the peak frequency deviation for 
STOs can be defined as 27rA/„ = 2z^rp|X„| for given 
harmonic frequency and modulation conditions. Due to 
the form of \Xn\^ it is expected to be linearly depen- 
dent on the modulation strength p and proportional to 
[(nwm)^ -I- (2rp)^]~^/^. The latter dependence is of par- 
ticular interest in terms of the modulation bandwidth 
(MBW) discussed in Section III. 

Several features of the spectrum can be readily iden- 
tified from the right-hand side of Eq. ([TU]). First we 
can see the shifted carrier frequency which, by virtue 
of the properties of convolution, shifts the base-band 
spectrum in frequency. Second, a nonlinear amplitude- 
modulation (NAM) spectrum arises from the power fluc- 
tuations. Third, a series of FM spectra whose harmonics 
expand with n can be seen to arise from the combined 
contributions of the phase modulation and the nonlin- 
early enhanced power fluctuations. In the case of a weakly 
nonlinear oscillator [v < 1), the FM scries of convolutions 
can be neglected, and the spectrum reduces to a NAM 
spectrum. For a particular applied fleld angle, the condi- 
tion J/ = is satisfied, and the spectrum reduces to pure 
AM4^ On the other hand, if the time-dependent term is 
neglected (C2 « 0), only the first harmonic coefficients 
of Eq. (jASp will be non-zero, and the PSD will take the 
form of pure FM4i 

The convoluting terms in Eq. (jlOp deflne the power 
of each harmonic. Although it is tedious to obtain a 
meaningful analytical expression, it can be inferred that 
the power of the harmonics is sensitive to the phase ipn- 
We can illustrate this by approximating each convolut- 
ing term to the second harmonic, and solving for the first 
upper and lower sideband power. By expanding the abso- 
lute values to first order, we obtain the power difference 
of the sideband A = 2Ji(/3i)(l + Ao)po, which is in gen- 
eral non-vanishing. This crude approximation does not 
represent the asymmetry of Eq. pH)) . which is affected 
up to the fifth harmonic, as shown in Appendix A. 

Summarizing this section, we have obtained a fairly 
complex spectrum by solving the set of Eq. ([S]). This 
complexity, along with the interdependence of several 
key parameters such as the modulation strength and fre- 
quency, limits the analytical insight that can be gained. 
However, the frequency shift of the carrier can be explic- 
itly obtained by considering a Fourier-series approach, 
and the sideband asymmetry arises as a consequence of 
the solution's harmonic-dependent phases. The solution 
of the spectrum leads to the definition of the harmonic- 
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FIG. 2. (Color online) (a) Frequency dispersion of the STO 
considered in the numerical calculations, as a function of the 
supercriticality ^. Beginning from ^ = 1, current-driven os- 
cillations are sustained. The frequency increases nonlinearly 
from the FMR frequency indicated by the horizontal line. 
Three operation points, OPl, 0P2, and OPS, are selected for 
the following results, (b-d) PSD as a function of the modula- 
tion index /3i and normalized frequency scale. The frequency 
shift is proportional to both the modulation strength and the 
curvature of (a), following Eq. ([9]) (solid white line over the 
carrier). 



dependent modulation index, from which the STO peak 
frequency deviation can be defined. 



III. NUMERICAL RESULTS 

In this section, we consider a first order expansion 
of the STOs parameters in order to evaluate the PSD 
of Eq. pUj) . The coefficients are calculated up to the 
tenth harmonic by solving the linear system of equations 
Eq. (|A3[) in matrix form. We consider the nanocontact 
geometry of Fig. [T] with the parameters ^oHapp = 1 T, 
/jqMs = 0.8 T, 7 = 28 GHz/T, a = 0.01, and = 100. 
Evaluating the STO parameters r+(p), r_(p), and uj{p) 
according to auto-oscillator theory;^ we obtain the fre- 
quency vs supercriticality dispersion shown in Fig. [5^. As 
expected, the frequency is equal to the FMR frequency 
at the threshold of the oscillation (^ ~ 1). 

Three operation points (OP), OPl, 0P2, and 0P3, 
were selected in Fig. [5^, such that their restoration 
rates were rp/27r = 11.2, 44.8, and 156.8 MHz, respec- 
tively. Each OP was modulated, and their PSDs, calcu- 
lated as the FFT of Eq. ([BT]), are shown in Fig. [IJb-d) 
as a function of /3i and a fixed modulation frequency. 
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FIG. 3. (Color online) (a-c) Carrier (black) and first upper 
(red) and lower (blue) sidebands for each OP as a function of 
the modulation index /3i. The sideband asymmetry is clearly 
visible, and is substantially higher at OPl due to the strong 
frequency shift. Panels (d-f) show the normalized power dif- 
ference of the sideband A as a function of both modulation 
strength and frequency. White lines denote the slice shown 
in (a-c). The colormaps show that the asymmetry increases 
together with both modulation parameters. 



/,„ = 100 MHz. The first-harmonic modulation index is 
used as a reference, since it is assumed to be the main 
contribution to the sideband power. In order to evalu- 
ate Eq. ()A3|) . the modulation strength is calculated from 
/?! , returning different ranges for each OP. The frequency 
axis is normalized to the modulation frequency and cen- 
tered on the carrier frequency. On this scale, all opera- 
tion points share the same number of visible sidebands, 
as well as the position of the carrier. It is observed that 
the frequency shift is more pronounced near the thresh- 
old (OPl). This is expected from Eq. (j9]), shown as the 
white lines in Fig. [^b-d) , since the modulation strength 
is higher near both the threshold and the curvature of 
Fig [2ja) . Such curvature-dependence has been observed 
in experiments^ away from the oscillation threshold, pre- 
sumably due to higher order nonlinearities. These ef- 
fects can be included in the present framework through 
Eq. (j6bp . as the derivative of the negative damping may 
be substantially different. A numerical fit of experimental 
data is not within the scope of the present paper, but we 
argue that such a fit could be achieved by measuring the 
auto-oscillator parameters, as has recently been shown in 
Rcf. 0. 

The power of the normalized carrier and of the first up- 
per and lower sidebands is shown in Fig. [31[a-c) for each 
OP and with modulation frequency fm = 100 MHz. The 
combined action of the frequency shift and the harmonic- 
dependent phase introduced by the Fourier series results 
in power asymmetry between the sidebands. This asym- 



5 



metry is notably higher for OPl in correlation with its 
similarly enhanced frequency shift. In the color plots of 
Fig. |3|d-f), the power difference of the first upper and 
lower sidebands A is shown as a function of /3i and /,„. 
The white horizontal line denotes the slice shown in the 
panels (a-c), respectively. Being a common feature, the 
sideband asymmetry increases together with both the 
modulation strength and the frequency. These dependen- 
cies can be qualitatively understood from Eq. pO)) and 
Eq. (|A3|) . A stronger fj, increases the power of the higher 
harmonic coefficients, so that they become non-negligible 
in the PSD. Thus, an increase in the asymmetry is ex- 
pected. On the other hand, the dependence on the modu- 
lation frequency arises from the relation between A„ and 
Bn- Roughly, one can approximate Ai = uim/'2TpBi, so 
that the complex variables in the PSD can be expressed 
as Xi = Bi{l + icL>m/2rp). From here, it is clear that as 
LJrn increases past the parameter 2Tp, the imaginary term 
becomes dominant or, in other words, the harmonic phase 
increases towards ±7r/2. Consequently, the power contri- 
bution of each harmonic to the sideband becomes heavily 
weighted by the phase, leading to enhanced asymmetry. 

Finally, we discuss the features of the peak frequency 
deviation A/„ in relation to the modulation bandwidth. 
The frequency dependence of A/ for each OP and fixed 
/i = 0.05 is shown in Fig. IH both by evaluating the first- 
harmonic modulation index, Eq. (|lip (solid lines), and by 
calculating the maximum of the instantaneous frequency 
from Eq. (|8b[) (circles) . The magnitude of the deviation is 
shown relative to the OP oscillator frequency, while the 
logarithmic scale is used to enhance the frequency de- 
pendence. Both approaches agree very well, so that the 
first harmonic modulation index suffices to perform the 
following calculations approximately. A practical conse- 
quence of this feature is that choosing a maximum mod- 
ulation strength does not guarantee the frequency ex- 
cursion estimated from the STO's frequency vs current 
characteristics. This is closely related to the modulation 
bandwidth of oscillators. 

The modulation bandwidth (MBW) gives a measure 
of the frequency range in which an oscillator has opti- 
mal modulation properties. A common criterion is the 
3 dB power attenuation that is likewise used to charac- 
terize filters. Since the power in a FM scheme depends 
on the Bessel functions, the methods used to measure 
the MBW rely on indirectly estimating the degradation 
of the linearity of the modulation index. In the present 
framework, the modulation index is defined from the cal- 
culation of the power spectrum (Appendix B), and one 
can directly estimate the MBW from Eq. pT|) . This can 
be accomplished by considering a vanishing modulation 
frequency and an arbitrary ^, so that the first-harmonic 
modulation index is The task is then to keep the ratio 
/i/a;„i constant while looking for /3/ = /3i/\^. Perform- 
ing this calculation gives the result that the MBW for 
STOs is 2Tp. Since the lack of linearity in our framework 
is given by 27rA/, one can now understand Fig. 2] as a 
STO transfer function under modulation having a char- 
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FIG. 4. (Color online) Peak frequency deviation A/ depen- 
dence on the modulation frequency. The solid lines have been 
approximately calculated from the first-order modulation in- 
dex, while the circles are have been calculated from the max- 
imum amplitude of the instantaneous frequency. The dashed 
lines represent the cut-off frequency for each OP while 
the dash-dotted black line is a guide for the eye, showing a 
0.1 %/decade slope. 



acteristic low-pass filter form with a 0.1 %/decade roll-off 
slope (dash-dotted black line) after the cut-off frequency 
(dashed vertical lines). 



IV. CONCLUSIONS 

The NFAM spectrum is obtained from the auto- 
oscillator general model by considering the power per- 
turbation created by the modulating signal. These power 
fiuctuations are enhanced by the STO's nonlinearity, re- 
sulting in a spectrum consisting of a series of convolu- 
tions between NAM and FM spectra. Its implicit depen- 
dence on the total damping parameter and nonlineari- 
ties suggests that the NFAM characteristics will exhibit 
sample-to-sample variation. In fact, such characteristics 
have been observed in experiments and their depen- 
dence on the modulation parameters qualitatively agrees 
with the results shown here. In contrast to the originally 
proposed model;^ we found a frequency dependence of 
the modulation index which leads to the estimation of the 
modulation bandwidth for STOs. Moreover, the STO un- 
der modulation has a low-pass filter behavior with a cut- 
off frequency given by 2rp, which usually depends on the 
operation point of a specific sample. Although the mod- 
ulation bandwidth for STOs has yet to be measured ex- 
perimentally, we expect a correlated degradation of STO 
modulation characteristics in frequency-dependent stud- 
ies. On the other hand, the Fourier series solution pro- 
posed in this paper can, in principle, reduce to a linear 
set of coupled equations any STO geometry described by 
the auto-oscillator general equation Eq. ([5]) perturbed by 
slow time- varying signal. 

Support from the Swedish Research Council (VR) 
is gratefully acknowledged. Johan Akerman is a Royal 
Swedish Academy of Sciences Research Fellow supported 
by a grant from the Knut and Alice Wallenberg Founda- 
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tion. 



Appendix A 

The Fourier series solution proposed in Eq. ([7]) is intro- 
duced into Eq. ()5a|) . The variable term can be expanded 
using trigonometric identities: 



cos{uj„it)Sp = Aq cos(cj,„t) 

A. 



(Al) 

sin(ti;„i(l + n)t) — sin(wm(l — n)t)] 



^-2 

n=l 

-;^[cos(i:j„(1 + n)t) + cos(a;„(l - n)t)] 



This series can be then rearranged by expanding the 
summation and identifying equal harmonics. Changing 
the index accordingly, we obtain 



COs{uJjnt)Sp = (Aq + y-) 



(A2) 



-[{2B2A0) cos(cj™<) - A2 sin(w™t)] 



11=2 
{Bn-1 



This expression treats the harmonics separately, so 
that Eq. ([5a|) can be solved by collecting harmonic terms. 
The resulting system of equations for the coefBcients is 





UJmAi ■ 





C2 



(A3a) 



li{Ci-C2Ao- ^B2)-2Tj,Bu (A3b) 



2r„Ai 



^2 A 



E 

n=2 



2TpBn + nCdrnA. 
C2 



Bn+l), 



(A3c) 
(A3d) 

(A3e) 



c 

+^l^{An-l- An+l). 

There are two possible solutions to the set of 
Eq. (|A3p . One can, for instance, express a system of 
equations depending only on the As, so that An (x 
An+2T An+i, An-i, An-2- This ncw system can easily be 
solved numerically using matrix algebra. The error will 
depend on the size of the matrix or, in other words, on 
the maximum harmonic where the series is truncated. 
The total error is shown in Fig. [SJa) for two frequencies 
of different orders of magnitude. The case for = 20 har- 
monics is used as a reference. It is observed that in both 
cases. A > 5 returns a total error smaller than 10~^ %. 
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FIG. 5. (Color online) (a) Total error as a function of the 
number of harmonics A'^ in the series, calculated from the lin- 
ear equations Eq. (jA3|l . For N > 5, the error is assumed to 
be negligible. Two modulation frequencies are shown as lim- 
iting cases. For intermediate frequencies, the error lies in the 
shaded area, (b-c) Error between the linear set of equations 
and the recursive method as a function of the modulation 
index The modulation frequencies for each case are 40 
and 400 MHz. The error introduced by the recursive method 
is generally small, but diverges close to the minimum of the 
sideband, which is indicated by vertical dashed lines. 



Another solution is grounded on the fact that the har- 
monic coefficients decay with n, so that we can assume 
in general that A„ oc A„_i, An-2- This approach can be 
implemented numerically as a recurrent series, and this 
is shown for different values of A^ and fm in Fig. [5{b-c). 
Here we observe that N = 5 also converges to a minimum 
error in this approximation. However, as /3i is swept out, 
it is clear the the error increases towards the minimum of 
the first sideband (indicated by a dashed line). It is note- 
worthy that the maximum error in Fig. Efc) is slightly 
shifted, which corresponds to the impact of higher har- 
monic terms. 



Appendix B 

The power spectral density (PSD) of the proposed so- 
lution is obtained from Eq. ([5]), which defines both time- 
dependent power and phase variations. The PSD is de- 
fined as 



PSD^Po* /(I + M * /(cos(0)) 



(Bl) 



where / is the Fourier transform. In the following we per- 
form each Fourier transform separately, and express the 
result as a convolution. The first term on the right hand 
side has the form of a NAM. The phase introduced by 
the sine function leads us to define the complex variable 
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Xn = Bn + iAn , SO that 

f{l + 6p)^{l + Ao)S{0) (B2) 
+ ^ -^S{nuj,n) + -^(5(-7ia;„i), 

n=l 

where the notation d{xQ) = 6{x — Xq) is used and (5 is the 
Kronecker delta. 

The second term on the right hand side of Eq. (|Bip is 
calculated by using Euler's formulae, and subsequently 
expanding in Taylor series. We obtain 

where a„ = 2;^rpyl„/(ri,a;„i), &„ = 2vTpBn/{nuJm) and 
c.c is the complex conjugate. The two summations can 
be multiplied by expanding and rearranging the terms. 
Defining 2x„ = 6„ + ia„, it becomes possible to rewrite 



the product terms as 

oo oo _ 

n(-)-EM)'w 

oo oo _ 

From this expression, we can identify the summa- 
tions on k as Bessel functions of order j and argument 
j3n = 1\/xnXn (upon rcuormalization of the coefficients 
by Pn/'^)- The product takes the form 

oo 

X{{-) = UPn) 
n=l 

+ E^^[(a^»er + (-Xne-i)^]. (B5) 

The Fourier transform of Eq. (|B5p is easily performed 
by expressing the result as a product of convolutions. By 
further evaluating Xn in terms of the complex variable 
Xn, and convoluting each harmonic contribution, we ob- 
tain Eq. (fTOj) . 
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